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Abstract. Two-dimensional generalization of the original peak finding algorithm suggested 
earlier is given. The ideology of the algorithm emerged from the well known quantum mechanical 
tunneling property which enables small bodies to penetrate through narrow potential barriers. We 
further merge this "quantum" ideology with the philosophy of Particle Swarm Optimization to get the 
global optimization algorithm which can be called Quantum Swarm Optimization. The functionality 
of the newborn algorithm is tested on some benchmark optimization problems. 
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1. Introduction. Some time ago we suggested a new algorithm for automatic 
photopeak location in gamma-ray spectra from semiconductor and scintillator detec- 
tors The algorithm was inspired by quantum mechanical property of small balls 
to penetrate through narrow barriers and find their way down to the potential wall 
bottom even in the case of irregular potential shape. 

In one dimensional case the idea was realized by means of finite Markov chain and 
its invariant distribution T" . States of this Markov chain correspond to channels of the 
original histogram. The only nonzero transition probabilities are those which connect 
a given state to its closest left and right neighbor states. Therefore the transition 
probability matrix for our Markov chain has the form 
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As for the transition probabilities, the following expressions were used 
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The number m is a parameter of the model which mimics the (inverse) mass of the 
quantum ball and therefore allows to govern its penetrating ability. 

The invariant distribution for the above described Markov chain can be given by 
a simple analytic formula |2] 
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where ui is defined from the normalization condition 
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= 1 . 



Local maximums in the original spectrum are translated into the very sharp peaks 
in the invariant distribution and therefore their location is facilitated. 

The algorithm proved helpful in uniformity studies of NaJ(Tl) crystals for the 
SND detector Another application of this "peak amplifier" , to refine the amplitude 
fit method in ATLAS Bg-mixing studies, was described in g]. In this paper we will 
try to extend the method also in the two-dimensional case. 

2. Two-dimensional generalization. The following two-dimensional general- 
ization seems straightforward. For two-dimensional n x n histograms the correspond- 
ing Markov chain states will also form a two-dimensional array Let u^-*^^ be 
a probability for the state to be occupied after fc-steps of the Markov process. 
Then 
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where Pij-im is a transition probability from the state to the state (l, m). We will 
assume that the only nonzero transition probabilities are those which connect a given 
state to its closest left, right, up or down neighbor states. Then the generalization of 
equations Hl.l|l and 11.2|l is almost obvious. Namely, for the transition probabilities 
we will take 
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We are interested in invariant distribution Uij for this Markov chain, such that 
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Unfortunately, unlike to the one-dimensional case, this invariant distribution can not 
be given by a simple analytic formula. But there is a way out: having at hand the 
transition probabilities Pij-irm we can simulate the corresponding Markov process 
starting with some initial distribution u^^\ Then after a sufficiently large number 
of iterations we will end with almost invariant distribution irrespective to the initial 



choice of u-. 
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For example, in the role of u^-^^ one can take the uniform distribution: 
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For practical realization of the algorithm, it is desirable to have precise meaning 
of words "sufficiently large number of iterations" . In our first tests the following 
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Fig. 2.1. The upper histograms represent the initial data in the contour and lego formats respec- 
tively. The middle histograms show the corresponding probability distribution after 258 iterations 
for the penetrating ability m = 3. The lower histograms represent the invariant distribution for the 
penetrating ability m = 30. 
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stopping criterion was used. One stops at fc-th iteration if the averaged relative 



difference between w^^-* and probability distributions is less than the desired 



accuracy e: 
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The performance of the algorithm is illustrated by Fig l2.1l for a 100 x 100 histogram 
representing three overlapping Gaussians with different widths. As expected, it works 
much like to its one-dimensional cousin: the invariant probability distribution shows 
sharp peaks at locations where the initial data has broad enough local maximums. 
Note that in this concrete example the one iteration variability e = 10^^ was reached 
after 258 iterations for m = 3 and after 113 iterations for m = 30. 

Convergence to the invariant distribution can be slow. In the example given 
above by Fig l2.1l the convergence is indeed slow for small penetrating abilities. If we 
continue iterations for m — 3 further, the side peaks will slowly decay in favor of the 
main peak corresponding to the global maximum. In the case of m = 3 it takes too 
much iterations to reach the invariant distribution. However, as Fig l2 . II indicates . the 
remarkable property to develop sharp peaks at locations of local maximums of the 

(k) 

initial histogram is already revealed by m-^- when number of iterations k is of the 
order of 300. 

One can make the algorithm to emphasize minimums, not maximums, by just 
reversing signs in the exponents: 
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This is illustrated by Fig l2.2l Here the initial histogram is generated by using a variant 
of the Griewank function 5 

.^ ^ - 100)2 + (y - 100)2 , ^^^^ J/-100 

(2.4) F{x,y)^^ ^-^^ L_cos(x-100)cosL_ + i. 

This function has the global minimum at a point x = 100, y = 100 and in the 
histogramed interval 50 < x < 150, 50 < y < 150 exhibits nearly thousand local min- 
imums. Many of them are still visible in the probability distribution for penetrating 
ability m = 3. But for m = 30 only one peak, corresponding to the global minimum, 
remains. 

3. Quantum swarm optimization. The discussion above was focused on two- 
dimensional histograms, while in practice more common problem is finding global 
optimums of nonlinear functions. The algorithm in the form discussed so far is not 
suitable for this latter problem. However it is possible to merge its ideology with the 
one of particle swarm optimization "S", T, H to get a workable tool. 

The particle swarm optimization was inspired by intriguing ability of bird flocks to 
find spots with food, even though birds in the flock had no previous knowledge of their 
location and appearance. "This algorithm belongs ideologically to that philosophical 
school that allows wisdom to emerge rather than trying to impose it, that emulates 
nature rather than trying to control it, and that seeks to make things simpler rather 
than more complex" 6 . This charming philosophy is indeed very attractive. So we 
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Fig. 2.2. The probability distribution for the Griewank function. Upper histograms for m = 3, 
lower histograms for m = 30. 

attempted to develop quantum swarm optimization - when each particle in the swarm 
mimics the quantum behavior. 

The algorithm that emerged goes as follows: 

• initialize a swarm of Up particles at random positions in the search space 

•^min — — •^max: ymin ^ U ^ Umax- 

• find a particle ii, in the swarm with the best position {xb,yb), such that the 
function F{x, y) under investigation has the most optimal value for the swarm 

at (xb.yb)- 

• for each particle in the swarm find the distance from the best position d = 

{x — Xf,)^ + (y ^ yb)"^ ■ For the best particle instead take the maximal value 
of these distances from the previous iteration (or for the first iteration take 

^ ~ \/ {p^max •^miri)'^ ~l~ iymax ymin)'^ ) ■ 

• generate a random number r uniformly distributed in the interval < r < 1 
and find the random step h = rd. 

• check for a bettor optimum the closest left, right, up and down neighbor states 
with the step h. If the result is positive, change % and {xb,yb) respectively. 
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Otherwise 

move the particle to left, right, up or down by the step h according to the 
corresponding probabilities of such jumps: 

PL = : ; : , PR 
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and 

Xu = min (x + h, Xmax): ~ max (a; /l, Xmin): 

(3.3) yu = min{y + h,ymax), yd = max (y - /i, y„i„). 

At last /g = 1, if optimization means to find the global maximum, and Ig = 
— 1, if the global minimum is searched. 

• do not stick at walls. If the particle is at the boundary of the search space, 
it jumps away from the wall with the probability equaled one (that is the 
probabilities of other three jumps are set to zero). 

• check whether the new position of the particle leads to the better optimum. 
If yes, change if, and {xb,yb) accordingly. 

• do not move the best particle if not profitable. 

• when all particles from the swarm make their jumps, the iteration is finished. 
Repeat it at a prescribed times or until some other stopping criteria are met. 

To test the algorithm performance, we tried it on some benchmark optimization 
test functions. For each test function and for each number of iterations one thou- 
sand independent numerical experiments were performed and the success rate of the 
algorithm was calculated. The criterion of success was the following 

r 10-3|a:„|, if \x,n\ > 10-3 

(o,^ |„ „ I < / 10-3 |yml, i/ |yml > 10-3 

^^■'^ '^^"^"^'-1 10-3, i/ |y™|< 10-3 ' 

where {xm, ym) is the true position of the global optimum and {xf,yf) is the position 
found by the algorithm. The results are given in the tablcl^ll. The test functions itself 
are defined in the appendix. Here we give only some comments about the algorithm 
performance. 



Table 3.1 

Success rate of the algorithm in percentages for various test functions and for various numbers 
of iterations. Swarm size Up = 20. 



Function 

AT 

Name 


Iterations 
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600 
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Chichinadze 


35.5 


97 


100 


100 


100 


100 


100 


100 


Schwefel 


99.4 


99.5 


99.8 


99.3 


99.2 


99.8 


100 


99.6 


Ackley 


100 


100 


100 


100 


100 


100 


100 


100 


Matyas 


88.9 


100 


100 


100 


100 


100 


100 


100 


Booth 


100 


100 


100 


100 


100 


100 


100 


100 


Easom 


93.6 


100 


100 


100 


100 


100 


100 


100 


LcvyS 


98.4 


99.5 


99.4 


99.3 


99 


99 


99.1 


99.5 


Goldstein-Price 


100 


100 


100 


100 


100 


100 


100 


100 


Griewank 


76.3 


99.7 


100 


100 


100 


100 


100 


100 


Rastrigin 


100 


100 


99.8 


99.9 


100 


99.9 


99.9 


100 


Rosenbrock 


43.6 


90.4 


99.8 


100 


100 


100 


100 


100 


Leon 


13.8 


52.1 


82 


91.6 


97.6 


99.1 


99.6 


99.8 


Giunta 


100 


100 


100 


100 


100 


100 


100 


100 


Beale 


99.7 


100 


100 


100 


100 


100 


100 


100 


Bukin2 


61.8 


84.4 


93.8 


97.8 


98.6 


99.3 


99.7 


99.8 


Bukin4 


99.6 


100 


100 


100 


100 


100 


100 


100 


BukinG 


0.2 


0.1 





0.2 





0.1 


0.2 


0.1 


Styblinski-Tang 


100 


100 


100 


100 


100 


100 


100 


100 


Zettl 


100 


100 


100 


100 


100 


100 


100 


100 


Three Hump Camel 


100 


100 


100 


100 


100 


100 


100 


100 


S chaff c^r 


8.2 


34.7 


60.7 


71.2 


77.8 


78.9 


80.4 


83.9 


LcxyVS 


lUU 


100 


lUO 


100 


100 


lOU 


100 


100 


McCormic 


100 


100 


100 


100 


100 


100 


100 


100 



For some test problems, such as Chichinadze, Ackley, Matyas, Booth, Easom, 
Goldstein-Price, Griewank, Giunta, Beale, Bukin4, Styblinski-Tang, Zettl, Levyl3, 
McCormic and Three Hump Camel Back, the algorithm is triumphant. 

Matyas problem seems easy, because the function is only quadratic. However it is 
very flat near the line x = y and this leads to problems for many global optimization 
algorithms. 

Easom function is a unimodal test function which is expected to be hard for 
any stochastic algorithms, because vicinity of its global minimum has a small area 
compared to the search space. Surprisingly our algorithm performs quite well for 
this function and one needs only about 100 iterations to find the needle of the global 
minimum in a haystack of the search space. 

Schwefel function is deceptive enough to cause search algorithms to converge in the 
wrong direction. This happens because the global minimum is geometrically distant 
from the next best local minima. In some small fraction of events our algorithm is also 
prone to converge in the wrong direction and in these cases the performance seems 
not to improve by further increasing the number of iterations. But the siiccess rate is 
quite high. Therefore in this case it is more sensible to have two or more independent 
tries of the algorithm with rather small number of iterations each. 

Rastrigin function is a multimodal test function which have plenty of hills and 
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valleys. Our algorithm performs even better for this function, but the success is not 
universal either. 

Rosenbrock function is on contrary unimodal. Its minimum is situated in a banana 
shaped valley with a flat bottom and is not easy to find. The algorithm needs more 
than 200 iterations to be successful in this case. Leon function is of the similar nature, 
with even more flat bottom and the convergence in this case is correspondingly more 
slow. 

Griewank, LevyS and LevylS are multimodal test functions. They are considered 
to be difficult for local optimizers because of the very rugged landscapes and very 
large number of local optima. For example, LevyS has 760 local minima in the search 
domain but only one global minimum and LevylS has 900 local minima. Test results 
reveal a small probability that our algorithm becomes stuck in one of the local minima 
for the Levy5 function. 

Giunta function simulates the effects of numerical noise by means of a high fre- 
quency, low amplitude sine wave, added to the main part of the function. The algo- 
rithm is successful for this function. 

Convergence of the algorithm is rather slow for Bukin2 function, and especially 
for the Schaffer function. This latter problem is hard because of the highly vari- 
able data surface features many circular local optima, and our algorithm becomes, 
unfortunately, often stuck in the optima nearest to the global one. 

At last, the algorithm fails completely for the Bukin6 function. This function has 
a long narrow valley which is readily identified by the algorithm. But the function 
values differ very small along the valley. Besides the surface is non-smooth in the valley 
with numerous pitfalls. This problem seems hopeless for any stochastic algorithm 
based heavily on random walks, because one has to chance upon a very vicinity of 
the global optimum to be successful. The non-stochastic component of our algorithm 
(calculation of jump probabilities to mimic the quantum tunneling) turns out to be 
of little use for this particular problem. 

4. Concluding remarks. The Quantum Swarm Optimization algorithm pre- 
sented above emerged while trying to generalize in the two-dimensional case a "quan- 
tum mechanical" algorithm for automatic location of photopeaks in the one dimen- 
sional histograms ^j. 

" Everything has been said before, but since nobody listens we have to keep 
going back and beginning all over again" [H]. After this investigation was almost 
finished, we discovered the paper ^HI by Xie, Zhang and Yang with the similar idea 
to use the simulation of particle- wave duality in optimization problems. However their 
realization of the idea is quite different. 

Even earlier, Levy and Montalvo used the tunneling method for global optimiza- 
tion but without referring to quantum behavior. Their method consisted in a 
transformation of the objective function, once a local minimum has been reached, 
which destroys this local minimum and allows to tunnel classically to another valley. 

We found also that the similar ideology to mimic Nature's quantum behavior in 
optimization problems emerged in quantum chemistry and led to such algorithms as 
quantum annealing and Quantum Path Minimization |1S| . 

Nevertheless, the Quantum Swarm Optimization is conceptually rather different 
from these developments. We hope it is simple and effective enough to find an eco- 
logical niche in a variety of global optimization algorithms. 

Appendix. Here we collect the test functions definitions, locations of their opti- 
mums and the boundaries of the search space. The majority of them was taken from 
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|14l 1151 , but we also provide the original reference when known. 
Chichinadze function |1K[I17| 

TT 1 / (v — 5)^^ 

F{x,y) = - 12x + 11 + 10 cos -a; + 8sin(57ra;) - exp \ ~— — 

-30 <x,y < 30, F™„(x, y) = i^(5.90133, 0.5) = -43.3159. 
Schweffel function ^21 

y) = -X sin - y sin i/jyl, 

-500 < X, y < 500, Frmn {x,y) = F(420.9687, 420.9687) = -837.9658. 
Ackley function U^l 

F[x, y) — 20[1 — e^°-2'\/0-5(a;^+y^)j _ g0.5[cos(27rx)+cos(27ra)] _|_ 

-35 < a:,j/< 35, F™„(a:, y) = ^^(0, 0) = 0. 
Matyas function ^B] 

F{x, y) = 0.26(x2 + y^) - 0.48a;?;, 
-10<x,y<10, F™„(x,2/) = F(0,0) = 0. 

Booth function 

y) = {x + 2y- if + [2x + y - hf 

-10<x,y<10, i^™„(a;,y) = F(l,3) = 0. 
Easom function [221 

y) = — cos a; cos y exp [— (x — tt)^ — (y — i")^], 

-100 < a;,y < 100, F™„(a;,y) = F(7r,^) = -1. 
Levy5 function jJS] 

5 5 

F{x, y) = cos [{i - l)a; + X! [(■?' + + ■?'] + 

i=i j=i 

+ {x + 1.42513)2 + {y + 0.80032)2, 

-100 < X, y < 100, F„i„(x, y) = F(-l. 30685, -1.424845) = -176.1375. 
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Goldstein-Price function 

F{x, y)= [\ + {x + y + 1)^(19 - 14a; + ix^ - Uy + &xy + 3^^)] x 



[30 + [2x - 3?/)2(18 - 32a; + 12a; 



-2<a;,y<2, F„„„(x, y) = F(0, -1) = 3. 
Griewank function JSj [T5] 

^^"^'^^ " ~200 cosa;cos--^ + 1, 

-100 < X, 2/ < 100, F™„(a;, y) = F(0, 0) = 0. 
Rastrigin function 

F{x, y) = a;2 + - 10 cos (27ra;) - 10 cos (27ry) + 20, 

-5.12 < X, y < 5.12, y) = F(0, 0) = 0. 

Rosenbrock function JS] 

i^(x,y) = 100(y-a;2)2 + (l-x)2, 



Leon function 



-1.2 < X, y < 1.2, F™„(x, y) = F(l, 1) - 0. 



F(x, y) = 100(y - x^f + (1 - a;)^ 



-1.2 < x,y < 1.2, F„„„(x,y) = F(l, 1) = 0. 

Giunta function 



^/ N • ('16 \ . 2 \ 1 . 

y ) = sm — x — 1 + sm — x — 1 H sm 

^ ''^^ * 15 / Vl5 y 50 



4(i^x-l 
15 



,16 \ , /16 \ 1 

' sm — y — 1 + sm — y — 1 H sm 

M5^ / V15 / 50 



, 16 

4(y^.-1 



0.6, 



-1 < X, y < 1, F„„(x, y) = i^(0.45834282, 0.45834282) = 0.0602472184 
Beale function 

F(x, y) = (1.5 - X + xyf + (2.25 - x + xy^ f + (2.625 - x + xy^)^ 



-4.5 < x,y < 4.5, ^^™„(x,y) = F(3,0) = 0. 
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Bukin2 function |H1 

F{x, y) = 100(y - O.Olx^ + 1) + 0.01(a; + 10)^ 

-15<x<-5, -3<zj<3, F„,^n{x,y)=F{~10,0) = 0. 
Bukin4 function [211 

F{x,y) = 100y2 + 0.01|x+10|, 

-15 < X < -5, -3 < y < 3, F^„(x,y) = F(-10,0) = 0. 
Bukin6 function 



F{x,y) = 100^/\y^^~0mx^ + 0.01\x + 10\, 

-15<x< -5, -3 < y < 3, F^„,{x, y) = F(-10, 1) = 0. 
Styblinski-Tang function 

F{x, y) = ^ [x'^ - 16a;2 + bx + IGy^ + 5y] , 

-5 < a;, y < 15, F™„(x, y) = i^(-2.903534, -2.903534) = -78.332. 
Zettl function ^ 

F{x, y) = (a;2 + _ 2xf + 0.25a;, 

-5 < a;, y < 5, F,„i„(a:, y) = F(-0.0299, 0) = -0.003791. 
Three Hump Camel back function Jl] 

x^ 

Fix, y) = 2x2 _ Q5^4 ^ _ ^ 2;y + y2, 

6 

-5<x,y<5, F™„(x,y) =F(0,0) = 0. 
SchafFer function [221 



F(x,y) = 0.5' «-V^^-0-5 



[1 + 0.001(x2 + y2)]2' 

-100 < X, y < 100, F™„(x, y) = ^^(0, 0) = 0. 
Levy 13 function ^1] 

F(x, y) = sin^ (37rx) + (x - 1)^ [l + sin^ (37ry)] + (y - 1)^ [l + sin^ (27ry)] , 

-10<x,y<10, F™„(x,y) = f^(l,l) = 0. 
McCormic function |27| 

F(x, y) = sin (x + y) + (x - y)^ - 1.5x + 2.5y + 1, 

-1.5 <x<4, -3<x<4 F™„(x,y) = F(-0.54719, -1.54719) = -1.9133. 
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